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Disclosed herein is a system for rapidly 
resolving position with centimeter— level accu- 
racy for a mobile or stationary receiver (4). 
This is achieved by estimating a set of pa- 
rameters that are related to the integer cycle 
ambiguities which arise in tracking the carrier 
phase of satellite downlinks (5, 6). In the pre- 
ferred embodiment, the technique involves a 
navigation receiver (4) simultaneously tracking 
transmissions (6) from Low Earth Orbit Satel- 
lites (LEQS) (2) together with transmissions (5) 
from GPS navigation satellites (1). The rapid 
change in the line-of-sight vectors from the re- 
ceiver (4) to the LEO signal sources (2), due 
to the orbital motion of the LEOS, enables the 
resolution with integrity of the integer cycle 
ambiguities of the GPS signals (5) as well as 
parameters related to the integer cycle ambigu- 
ity on the LEOS signals (6). These parame- 
ters, once identified, enable real-time centime- 
ter-level positioning of the receiver (4). 
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A System Using LEO Satellites for 
Centimeter-Level Navigation 



This invention was reduced to practice with support from NASA under contract number 
NAS8-39225. The U.S. Government has certain rights in the invention. 

I. Cross-Reference to Related Applications 

This application claims priority from U.S. Provisional Patent Application 60/041,184 filed 
03/21/97, which is incorporated herein by reference. 

II. Background 

Conventional satellite-based positioning techniques are based on the use of special navi- 
gation signals transmitted from several navigational satellites. In the global positioning 
system (GPS), for example, a constellation of GPS satellites transmit LI and L2 carrier 
signals modulated with C/A and P code signals. By measuring these code -signals, a user 
receiver can determine its position to an accuracy of several meters. 

To determine the user position with higher accuracy, a differential technique can be used. 
A reference receiver having a known position also measures the code signals and calculates 
its position. The reference receiver then calculates a differential correction by comparing 
its known position with this calculated position, and transmits this correction to the user 
receiver. Assuming the user receiver is near the reference station, it can use the differential 
correction data to improve the accuracy of its position estimate down to approximately 
1 meter. 

Various proposed techniques provide positioning accuracy on the order of 1 cm. In addi- 
tion to measuring the code signals from the GPS satellites, these techniques use carrier 
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phase measurements of the signals from the GPS navigational satellites. Typically, this 
carrier phase positioning technique uses differential carrier phase correction data from a 
reference station in order to improve performance. There is a significant difficulty inher- 
ent to this technique, however. When tracking a carrier signal of a navigational satellite 
transmission, one is able to directly measure the phase of the signal, but one cannot de- 
termine by direct measurement how many complete integer cycles have elapsed between 
the times of signal emission and reception. The measured carrier signal thus has an in- 
herent integer cycle ambiguity which must be resolved in order to use the carrier phase 
measurements for positioning. Consequently, much research in the art of satellite-based 
positioning has focused on resolving these cycle ambiguities in carrier phase measurements 
of GPS satellite signals. 

MacDoran and Spitzmesser [1] describe a method for deriving pseudoranges to GPS satel- 
lites by successively resolving integers for higher and higher signal frequencies with mea- 
surements independent of the integers being resolved. The first measurement resolves the 
number of C/A code cycles using a Doppler range; these integers provide for independent 
measurements to resolve the number of P code cycles, and so on for the L2 and Ll carriers. 
This technique, however, assumes exact correlation between satellite and user frequency 
standards (i.e., the user requires an atomic clock), and provides no means of correcting 
for atmospheric distortions. 

A similar technique, called dual- frequency wide-laning, involves multiplying and filtering 
the L2 and Ll signals from a GPS satellite to form a beat signal of nominal wavelength 
86 cm, which is longer than either that of the Ll signal (19 cm) or the L2 signal (24 
cm). Integer ambiguities are then resolved on this longer wavelength signal. Since the 
L2 component is broadcast with encryption modulation, however, this technique requires 
methods of cross- correlation, squaring, or partially resolving the encryption. These tech- 
niques are difficult to implement and have low integrity. 

Hatch [2] describes a technique for resolving integer ambiguities using measurements from 
redundant GPS satellites. Initial carrier-phase data is collected from the minimum num- 
ber of GPS satellites needed to resolve the relative position between two antennae. From 
these measurements, a set of all possible integer combinations is derived. Using carrier 
phase measurements from an additional GPS satellite, the unlikely integer combinations 
are systematically eliminated. This technique is suited to the context of attitude deter- 
mination where both receivers use the same frequency standard and the distance between 
the antennae is fixed. This approach, however, is ill-adapted for positioning over large dis- 
placements, where the initial set of satellites is four and the distance between the receivers 
is not known a priori - the technique is then extremely susceptible to noise, and compu- 
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tationally intensive. Knight [3] details an approach similar to that of Hatch, except that 
a more efficient technique is derived for eliminating unlikely integer combinations from 
the feasible set. Knight's technique also assumes that the two receivers are on the same 
clock standard. 

Counselman [4] discloses a technique for GPS positioning that does not resolve integer 
cycle ambiguity resolution but rather finds the baseline vector between two fixed antennae 
by searching the space of possible baseline vectors. The antennae track the GPS satellite 
signals for a period of roughly 30 minutes. The baseline is selected that best accounts 
for the phase changes observed with the motion of the GPS satellites. This technique, 
however, assumes that the baseline vector remains constant over the course of all the 
measurements during the 30 minute interval, and is therefore only suitable for surveying 
applications. Moreover, it also assumes that the clock offset between user and reference 
receivers remains constant over the 30 minute measurement interval. 

A motion- based method for aircraft attitude determination has been disclosed by Cohen 
[5]. This method involves placing antennae on the aircraft wings and tail, as well as a 
reference antenna on the fuselage. The integer ambiguities between the antennae can be 
rapidly resolved as the changes in aircraft attitude alter the antenna geometry relative to 
the GPS satellite locations. This approach, however, is limited to attitude determination 
and is not suitable for precise absolute positioning of the aircraft itself. 

Current state-of-the art kinematic carrier phase GPS navigation systems for absolute 
positioning have been disclosed by Cohen [6, 7, 8, 9] and by Pervan [10, 11]. These 
systems achieve rapid resolution of cycle ambiguities using ground-based navigational 
pseudo satellites (pseudolites) which transmit either an additional ranging signal (Doppler 
Marker) or a signal in phase with one of the satellites (Synchrolites). .Although this 
approach rapidly achieves high precision absolute positioning, it provides high precision 
and integrity only when the user moves near the ground-based pseudolites. In addition, 
the pseudolites are expensive to maintain. 

Therefore, each of the existing techniques for satellite- based navigation suffers from one 
or more of the following drawbacks: (a) it does not provide centimeter- level accuracy, (b) 
it does not quickly resolve integer cycle ambiguities, (c) it is not suitable for kinematic 
applications, (d) it provides only attitude information and does not provide absolute 
position information, (e) it does not have high integrity, (f ) it requires the deployment and 
maintenance of pseudolites, (g) its performance is limited to users in a small geographical 
area near pseudolites. or (h) it requires the user receiver and/or the reference receiver to 
have an expensive highly stable oscillator. 
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III. Summary of the Invention 

In view of the above, it is an object of the present invention to provide a system and 
method for centimeter-level kinematic positioning with rapid acquisition times and high 
integrity. In addition, it is an object of the invention to provide such a method that 
does not depend on additional signals transmitted from nearby navigational pseudolite 
transmitters, and does not require a highly stable oscillator such as an atomic clock. 
It is further an object of the invention to provide a navigation system requiring only 
carrier phase information. These together with other objects and advantages will become 
apparent in the following description. 

In order to obtain high- integrity estimation of integer cycle ambiguities, carrier-phase mea- 
surements must be made for a time interval long enough that the displacement vectors 
between the user and the signal sources undergo substantial geometric change. Surpris- 
ingly, the present inventors have discovered a method and system for fast acquisition, 
high integrity, kinematic carrier-phase positioning using non-navigational signals from 
low earth orbit (LEO) satellites which are not necessarily intended for navigational use. 
The short orbital periods of these LEO satellites provide the required change in geometry 
for resolution of cycle ambiguities with high reliability in a few minutes. The technique, 
therefore, provides fast acquisition, high precision and high integrity without depending 
on signals from ground-based pseudolites in close proximity to the user. In addition, the 
technique has the advantage that it does not require that the LEO satellites have any 
special features (e.g. atomic clocks or the ability to transmit navigational signals). 

Remarkably, the present inventors have discovered a system and method for satellite- based 
navigation using signals from non- navigational satellites in low earth orbit. Beginning only 
with an estimate of the user clock offset, high precision kinematic positioning is provided 
using only carrier-phase signals transmitted from earth orbiting satellites. By using signals 
from at least one LEO satellite, high integrity and fast acquisition is provided. The other 
signals may be from other satellites, including high earth orbit navigational satellites, 
or from any other space-based or earth-based sources, including pseudolites and other 
earth-based transmitters. Only the carrier phase signals from these other sources are 
required. 

In a preferred embodiment, an initial estimate of the user position and clock offset is 
provided by conventional code-phase differential GPS techniques. In addition, differential 
carrier phase measurements are used in order to eliminate errors caused by atmospheric 
phase distortion, satellite ephemeris deviations, and deliberate corruption of the satellite 
signals. As will become apparent, however, the fundamental technique of using non- 
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navigational carrier signals from LEO satellites for resolving integer cycle ambiguities in 
a navigational system is not limited to these specific implementations. In alternative 
embodiments, for example, an initial clock offset may be estimated by any combination of 
known techniques for navigation, including anything from sophisticated earth-based radio 
navigation to simply calibrating the user receiver to a known reference. 

In the preferred embodiment, centimeter-level positioning is provided by combining the 
navigational data available from GPS satellites with the non-navigational carrier phase 
data available from LEO satellites. In addition, the non-navigational carrier phase data 
from the GPS satellites is used. The method is robust to frequency- dependent phase- 
lags in the navigation receivers, as well as to instabilities in the crystal oscillators of the 
satellites, and of the receivers. We describe how the general technique can be applied to 
a variety of different satellite communication configurations. Such configurations include 
satellite transmission of multiple beams with different phase-paths, and bent-pipe archi- 
tectures where the uplink signal is frequency- converted by the satellite and retransmitted. 

Generally, in one aspect of the invention a user device is provided for satellite- based nav- 
igation. The device comprises at least one antenna for coupling to signals transmitted 
from a set of satellites. The set of satellites includes a set of LEO satellites that do not 
necessarily transmit navigational information. A receiver in the device tracks the signals 
to obtain carrier phase information comprising geometrically diverse carrier phase infor- 
mation from the LEO satellites. A microprocessor in the device calculates the precise 
position of the user device based on the carrier phase information and an initial estimate 
of the device clock offset. In a preferred embodiment, the device calculates an initial esti- 
mate of position and clock offset from code phase information derived from navigational 
signals transmitted by navigational satellites. In addition, the preferred embodiment uses 
reference carrier phase information transmitted from a reference station to improve the 
accuracy of the position estimate. 

In another aspect of the invention, a satellite- based navigation system is provided. The 
system comprises a set of satellites, including LEO satellites, that transmit carrier signals, 
a reference station, and a user device. The reference station samples the carrier signals to 
obtain reference carrier phase information which is then transmitted to the user device over 
a communication link. In addition to receiving the reference carrier phase information, 
the user device directly tracks the carrier signals to obtain user carrier phase information 
from the set of LEO satellites. The user device then calculates its precise position based 
on the reference carrier phase information and the user carrier phase information. The 
calculation uses the geometrically diverse reference and user carrier information from the 
LEO satellites to quickly resolve parameters related to the integer cycle ambiguities in the 
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reference and user carrier phase information. In a preferred embodiment, the calculation 
of the user position is based on an initial estimated clock offset and position calculated 
from navigational code phase signals transmitted from a set of navigational satellites. 
Preferably, the reference station also transmits differential code phase correction data to 
the user to improve the accuracy of the initial estimate. 

In another aspect of the invention, a method is provided for estimating a precise position 
of a user device in a satellite- based navigation system. The method comprises transmit- 
ting carrier signals from a set of satellites > wherein the set of satellites includes a set of 
LEO satellites; tracking at a user device the carrier signals to obtain user carrier phase 
information comprising geometrically diverse user carrier phase information from the set 
of LEO satellites; and calculating the precise position of the user device based on an 
initial position estimate and the user carrier phase information, wherein the calculation 
uses the geometrically diverse user carrier information from the set of LEO satellites to 
quickly resolve integer cycle ambiguities in the user carrier phase information. In a pre- 
ferred embodiment, the method includes tracking at a reference station the carrier signals 
to obtain reference carrier phase information comprising geometrically diverse reference 
carrier phase information from the set of LEO satellites. The reference carrier phase in- 
formation is then used to improve the accuracy of the position calculation. In a preferred 
embodiment, the method further comprises estimating an approximate user position and 
clock offset using code phase signals received from a set of navigational satellites. Prefer- 
ably, differential code phase techniques are used to improve the accuracy of the initial 
estimate. The preferred embodiment of the method also includes additional advanta- 
geous techniques such as: compensating for frequency dependent phase delay differences 
between carrier signals in user and reference receiver circuits, reading navigation carrier 
information and LEO carrier information within a predetermined time interval selected in 
dependence upon an expected motion of the user receiver and the LEO signal sources, cal- 
ibrating LEO oscillator instabilities using navigation satellite information, compensating 
for phase disturbances resulting from a bent pipe LEO communication architecture, com- 
pensating for oscillator instabilities in the user and reference receivers, predicting present 
reference carrier phase information based on past reference carrier phase information, and 
monitoring the integrity of the position calculation. 

* 

IV. Brief Description of the Drawing Figures 

Fig. 1 shows an operational overview of a preferred embodiment of the invention. 

Fig. 2 shows some possible methods of conveying the LEO ephemerides to the user in a 
system according to the invention. 
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Fig. 3 shows a receiver architectural overview according to a preferred embodiment of the 
invention. 

Fig's 4a and 4b show two different mixing, filtering and sampling schemes according to a 
preferred embodiment of the invention. 

Fig. 5 shows tracking and phase counting assemblies for the mixing scheme of Fig. 4a, 
according to a preferred embodiment of the invention. 

Fig. 6 shows tracking and phase counting assemblies for the mixing scheme of Fig. 4b, 
according to a preferred embodiment of the invention. 

Fig. 7 shows tracking and phase counting assemblies for the mixing scheme of Fig. 4b, 
and a phase latch architecture, according to a preferred embodiment of the invention. 

Fig. 9 shows a microprocessor block diagram for a user receiver according to a preferred 
embodiment of the invention. 

Fig. 10 shows a microprocessor block diagram for a reference receiver according to a 
preferred embodiment of the invention. 

Fig. 11 is a conceptual illustration of lattice points defined by a lattice basis G in accor- 
dance with a preferred embodiment of the invention. 

Fig, 12 shows the beam arrangement for the Globalstar S-band downlink. 

Fig. 13 shows the operational overview of a method using Globalstar satellites according 
to a preferred embodiment of the invention. 

Fig. 14 is a graph of the fractional availability of the Globalstar Constellation above 10° 
elevation. 

Fig. 15 is a graph of the availability of RAIM for protection radii of 110 cm, assuming 
GPS alone and GPS augmented with Globalstar, according to a preferred embodiment of 
the invention. 

Fig. 16 shows a reference receiver and transmitter architectural overview in accordance 
with a preferred embodiment of the invention. 

Fig. S is a block diagram of one channel of a tracking module designed for the CDMA 
signal of equ. (1), according to a preferred embodiment of the invention. 

Fig. 17 is a graph of average standard deviations in radial position errors for a mobile 
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user, at 5km and I km from the reference station, according to a preferred embodiment 
of the invention. 

Fig. 18 is a graph of the evolution of the Probability of Selecting the Correct Integer Set, 
in accordance with a preferred embodiment of the invention. 

V. Detailed Description 

For the purposes of this description, we assume that an initial estimate of position and 
clock offset is derived from satellite navigational signals and that the navigation satellites 
employed are the Navstar satellites of the GPS constellation. It will be appreciated that 
other present or future navigation satellites may be used in other embodiments of the 
invention and other navigation techniques may be used to estimate an initial position and 
clock offset. We assume also the use of differential techniques involving a single reference 
receiver and a single user receiver which requires real-time position information. It will be 
apparent, however, that these differential techniques are not necessary in order to practice 
the invention. 

Fig. 1 presents a schematic overview of the system. The central components of the system 
are the Navstar GPS satellites la-d, the LEO satellites 2a, b, the user receiver 4 and the 
reference receiver 3. The user 4 and reference 3 receivers track the absolute carrier phase 
of the Navstar satellite signals 5a-d together with the absolute carrier phase of multiple 
LEO satellite signals 6a,b. By the term absolute, we mean that the phase measurement 
is accumulated over time and is not modulus 2n. The motion of the LEOS 2a,b causes 
rapid change in the angles - shown in the figure as B\ and 6 2 - between the baseline 7 
from user antenna 17a to reference antenna 17d, and the line of site vectors from user 
antenna 17a to satellites 2a,b. The rapid change in the line-of-sight vectors enables the 
user receiver 4 to resolve the integer cycle ambiguities on the Navstar satellite signals 5 a- 
d as well as parameters related to the integer cycle ambiguities on the LEO signals 6a,b, 
and consequently to position itself with cm-level precision with respect to the reference 
receiver 3. 

We describe below the steps involved in the general technique. It would be clear to one 
skilled in the art how the order of these steps would change, or how some steps would be 
adjusted for a static user, or an attitude determination problem where user and reference 
receivers are driven by a common oscillator and are separated by a constant distance. 

•The reference 3 and user 4 receiver obtain up-to-date satellite ephemeris information. 
•Both reference 3 and user 4 receiver measure the code phase delay on the signals 5a- 
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d transmitted by the GPS satellites. This measurement is known in the art as raw 
pseudorange. 

•Based on the code phase measurement, the user 4 and reference 3 correlate their clocks 
to within 1 fisec of GPS time (see Parkinson [12]). 

•In a preferred embodiment, the reference receiver 3 calculates differential corrections for 
the code-phase measurements and conveys them to the user receiver 4 over a communi- 
cation link 8. The user receiver 4 then positions itself with meter-level accuracy relative 
to the reference receiver 3 using the differentially corrected code-phase measurements. 

•The user 4 and reference 3 receiver simultaneously track the absolute carrier phase of 
GPS satellite signals 5a-d and LEO satellite signals 6a,b over an interval of time. 

•If necessary, the reference receiver calibrates the LEO satellite oscillators using the tech- 
nique described in section (VII. C). 

•The reference receiver 3 conveys to the user its carrier phase measurements and measure- 
ment correction data over the communication link 8. 

•The user corrects for deterministic errors in the carrier phase measurement. 

•The user employs the data- reduction technique described in section (VLB) to identify the 

integer cycle ambiguities on the Navstar satellites signals 5a-d, and parameters related 

to the integer cycle ambiguities on the LEOS signals 6a,b. 

• Once these parameters are identified, the user receiver 4 is able to position itself in real- 
time, with centimeter-level accuracy relative to the reference receiver 3. 

V.A. Data Links in the System which have Many Possible Implementations 
V.A.I Communicating the LEO Satellite Ephemerides 

The navigation capability of the system hinges on the user knowing the location of the 
LEOS 2a,b to reasonably high accuracy - see section (VIILD). Knowledge-of the satellite 
position as a function of time can be obtained via the satellite ephemeris data. This 
ephemeris data consists of several parameters for each satellite, describing the satellite 
orbit and changes in that orbit over time. To maintain the desired levels of accuracy, the 
ephemeris data must be updated by the user roughly daily. Fig. 2 displays the different 
mechanisms by which the ephemeris data may be conveyed to the user. Any sequence 
of arrows leading from the satellite operations and control center - SOCC 9, or tracking 
station 10 to the user 4a-c is possible. The position sensing for the ephemeris data can 
be achieved either by position sensors on the satellites 2a-c, such as GPS receivers, or by 
a tracking station 10 processing Doppler information from ground receivers at surveyed 
locations to calculate the LEOS' 2a-c orbital parameters. Whether the information is 
attained from the SOCC 9, or from a separate tracking station 10, the data is conveyed 
to an ephemeris data provider 11 who must make the information accessible to the user 
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4a-c. A simple implementation connects the reference 3 to the ephemeris provider 11 via 
line modems 12a,c over a regular land telephone line 12b. Alternatively, the information 
can be obtained by the reference over a LEO satellite data link 13a-b. The reference 
would then convey this information to the user 4a-c via a data Link from Reference to 
User - LRU 8. Another embodiment has the LEO satellite 2b broadcasting ephemeris data 
14a-e on a dedicated broadcast channel 14b-e which is received by both the reference 
3 and the users 4a-c. The techniques for deriving ephemeris data from satellite position 
sensors, or from Doppler-tracking satellites, are well understood in the art, as are the 
implementations of these different methods for conveying that information to the user. 

V.A.2 The data Link from Reference to User (LRU) 

For a mobile user, the LRU 8 is implemented with a real-time radio connection 15a-d. 
For an attitude determination problem, the LRU 8 could be implemented with a real- 
time cable connection, and for static surveying applications, the LRU could be off-line. 
For the mobile user the LRU 8 will be implemented using the same basic frequency and 
modulation scheme of either the GPS signals 5a-d, or the LEO satellite signals 6a,b. In 
this way, existing receiver front-end hardware - see Fig. 3 - can be used for the LRU 8. 
The radio LRU could be implemented by transmitting a signal 15a,b directly to the user 
4a, b via a reference station transmitter 16, or by using an existing LEO satellite data 
link 15c, d. 

The central function of the LRU 8 is to convey carrier phase measurements made at 
the reference station 3 to the user 4. If the user receiver 4 knows the location of the 
reference station 3, it is able to accurately predict the phase measurements that the 
reference receiver 3 will make- for up to a few seconds. Consequently, the LRU 8 needs 
to be active every few seconds over the duration of navigation. In addition to the carrier 
phase information, the LRU 8 can convey to the user 4: 

•The satellite ephemerides. This would be necessary for an implementation where the 
ephemerides are not known to the user, and cannot be obtained by a directly broadcast 
satellite data link 14a-e. 

•An estimate of the reference receiver's clock offset. This can be used to correct the 
differential measurement as described in section (VI. A). 

•Differential corrections for code-phase measurement. In order to achieve an initial differ- 
ential position estimate which is accurate to within meters, the reference 3 would send 
to the user 4 a set of corrections for the range measurements to improve the code phase 
performance of GPS. The technique of deriving these corrections is well understood in 
the art (see, e.g., Parkinson [13]). 
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•The position of the reference station antenna 17d. The data reduction technique depends 
upon the user knowing the rough location of the reference station antenna. The user 
receiver 4 finds the centimeter-level position between its antenna 17a-c and the reference 
receiver antenna 17d. Therefore, the universal accuracy of the user's derived position 
depends on the reference antenna 17d position information. 

•Status information on the satellites la-d,2a-c. The reference station can also send the 
user information about the health and signal characteristics of each of the satellites being 
tracked. 

•Error correction information. This information is employed by the user to minimize the 
residual errors of differential phase measurement, due for example to ionospheric and 
tropospheric delays, and severe satellite oscillator instabilities. 

V.B. Description of a combined GPS and LEO receiver 

Fig. 3 illustrates the essential components for a reference 3 or a user 4 receiver, assuming 
that transmissions are receivable from the Navstar Constellation (N) as well as LEO 
constellations (jL 1? L2, £3). Only one LEO constellation is necessary for the invention, 
however 3 LEO constellations are assumed for this diagram in order to illustrate the 
scaleable nature of the invention. While the overall structure of this receiver cannot be 
substantially altered, it should be recognized that the individual modules identified in 
Fig. 3 can be implemented in a variety of ways. 

The antennae subsystem 17 must be sensitive to the downlink bands of each constellation 
being tracked. We assume in our description of the data reduction technique that the an- 
tenna subsystem consists of a single antenna which is resonant at the relevant frequencies 
and has phase centers 21a-d for the respective bands which are separated by only a few 
millimeters, (see e.g. Long [14] and Zhong [15]). This assumption simplifies the descrip- 
tion, but is not essential to the invention. If an antenna subsystem is used with phase 
centers which are substantially separated for the different frequencies, the relative posi- 
tion of the phase centers would simply be used to correct the carrier phase measurements 
- equ. (13) - to take this separation into account. This would involve either assuming an 
orientation of the antenna subsystem 17, or for attitude- determination problems, itera- 
tively making the correction based on the estimated orientation of the vehicle on which 
the antenna subsystem is mounted. 

The antenna feeds a low noise amplifier (LNA) 22 which has gain over the full bandwidth 
of the signals being tracked. Of course, separate LNA's could be used for each of the con- 
stellations being tracked, namely on the iV, LI, L2 and LZ paths 18, 19a-c respectively. 
We will examine the receiver path for one of the constellations, Ll 19a, since all paths 
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are similar at this level of description. The signal from the LNA 22 is bandpass filtered 
with an rf bandpass filter 23 and then downconverted by mixing with a locally generated 
rf frequency / r/L1 24. An if bandpass filter 25b at the mixer 26b output removes the 
unwanted upper sideband. A microprocessor-controlled automatic gain control 29 then 
adjusts the magnitude of the signal to achieve optimal use of the available sampling bits. 
For a situation where a signal is hidden in bandpass noise of substantially larger power - 
such as arises with spread spectrum modulation methods - the SNR forfeited by sampling 
with one bit is roughly 1.96 dB. Therefore, it is possible to implement this system with 
only one bit, in which case the controllable gain of the AGC is unnecessary. Another 
potential variation on this architecture is to place a low noise AGC at the input to the rf 
bandpass filter 23,30-32 for each of the receiver paths, to substitute for the LNA 22 and 
AGC's 29,33-35. 

After amplification by the AGC 29, the signal 53 enters an if mixing stage 27, followed 
by a filtering stage 36 and a sampling stage 37. The signal is mixed in an if mixer 27 with 
a locally generated if frequency /,/ L1 28. The exact means of downconverting, filtering 
and sampling the signal at the if stage varies substantially with the signal structure and 
designer's preference. Fig's 4a and 4b illustrate two possible means of implementing the 
if mixing, filtering and sampling. Both of these schemes would be possible for Binary- 
Phase- Shift- Keyed (BPSK) and Quadrature-Phase-Shift-Keyed (QPSK) signals. In the 
scheme of Fig. 4a, only a single mixer 40 and filter 41 are used, and the sampling section 
42 outputs a single sample 43. If the incoming signal 53 had quadrature modulation, 
Fig. 5 describes a possible architecture for a tracking assembly 44b preceded by the 
scheme of Fig. 4a. The tracking modules 48a-c could mix the incoming signal 43 with 
both an in-phase and quadrature component of the output of the Numerically Controlled 
Oscillators 49a-c, in order to isolate the in-phase and quadrature modulation. In the 
scheme of Fig. 4b, both an in-phase 45a and quadrature 45b mixer are used, each of 
which outputs to a separate filter 46a, b and sampler 47a,b. If the incoming signal 53 has 
quadrature modulation, Fig. 6 describes a possible architecture for a tracking assembly 
44b preceded by the scheme of Fig. 4b. The tracking modules 52a-c of the tracking 
assembly could mix the incoming I 50a and Q 50b signals with a single in-phase output 
of the Numerically Controlled Oscillator 51a-c to isolate the in-phase and quadrature 
modulation. 

The bandwidth of the filters 41, 46a,b is chosen to accommodate the full bandwidth of 
the signal being tracked, which is shifted by the offset frequency f 0 = f Lx — f r j Ll — fij Ll - 
These filters 41, 46a,b could be lowpass or bandpass, depending on / Q . The sampling rate 
should be roughly 5 to 10 times the highest frequency component in the signal and the 
number of bits sampled could vary from 1 to 16 bits, depending on the signal structure, 
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the SNR, the desired robustness to interference, and hardware costs. 

Since the tracking assemblies of Fig. 5 and 6 have essentially similar structure, we will 
consider in more detail only the architecture displayed in Fig. 6. We assume the if 
mixing is performed using the scheme of Fig. 4b. The thick black lines represent digitized 
I 50a and Q 50b samples. The incoming samples are latched 72a, 72b and input 
to S tracking modules, where 5 is the maximum number of signals from a particular 
satellite constellation that one seeks to track. Each tracking module 52a-c tracks one 
satellite downlink signal by means of a phase-locked loop. Many different techniques 
for implementing the tracking modules 52a-c are known in the art. Whichever tracking 
module architecture is employed, it involves an oscillator which is phase-locked to the 
phase of the incident signal. The preferred embodiment employs a numerically controlled 
oscillator 51a-c in the phase-locked loop. In order to analyze phase- tracking behavior, we 
present one possible design of a tracking module for a generic CDMA signal of the form 

s(t) = -^D(z)C/(Ocos(u?t + <£(x)) 
v 2 

+ -^=D(t)C Q {t)sin(u>t + <t>(i)) + n{t) (1) 
v2 

where D{t) refers to the outer data sequence modulated on both the in- phase and quadra- 
ture signals. Ci(t) and Cq(0 are respectively the in-phase and quadrature spreading 
sequences. n(t) Represents thermal input noise, which is assumed to be normally dis- 
tributed, of zero mean and of uniform spectral density N 0 . A tracking module designed 
for the signal structure of equ. (1) is illustrated in Fig. 8. Ignoring the effect of front-end 
gains equally applied to signal and noise, the in-phase 98a and quadrature 98b digital 
signal entering the tracking module can be described: 

h = ^y=D k C Jk cos(<l> k ) + I n k 
A 

^y=DkC Q kCOs(<t>k) + Qnk 

E{P nk } = = ^ (2) 



Qk = 



where B c is the pre-correlation signal bandwidth, determined by the filters 46a, b. The 
upper sideband emerging from the carrier mixers 100a, b has frequency ~ 2/ 0 and is 
rejected by the accumulator lOla-d which has effective bandwidth ^, where T is the 
period of the inner codes Cj(£) and Cq(1). Consequently, we consider only the lower 
sideband of the mixer outputs 85a, b: 

A ^ 

Ilk = ^y=D k Cl k COs(<f> k — <t>rk) + link 
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A 

Qlk — -^D k CQhCOs((j>k ~ (Prk) + Qlnk 

B{tLl = E{Q\ nk } « ^ (3) 



Each of these signals is then mixed with a prompt inner code 102a,b and a tracking 
inner code 103a,b, which consists of the difference between early and late code replicas, 
separated by some number of chips, <f, where d < 2. The prompt accumulator 101b, c 
outputs can be described: 

A M 

—^R(Ti)Di ]T COs(<j> k - (prk) + hni 

A ™ 

-j=R(Ti)Di 2J Sin(<f> k - <p rk ) + Q2ni 

v2 fc=1 



Q 2 . = 



* {jl} - * {<?y * ^ (4) 

where -R(t,-) is the cross-correlation between the incoming code and the generated code 
for a time misalignment of r t . R(t) 1 once code- lock is achieved. For an inner code 
period, T, and sample rate / 5 , the summation is typically over M = T * f s samples. This 
number may be varied to accommodate code Doppler. Assume that over one correlation 
period, T, the NCO 105 suffers a constant frequency error, A/, <C ^, and a constant 
phase error, A<f>{. We can then treat the summation in equ. (4) as a continuous integral, 
and find 

A 

E {I 2 i} « -j^Disinc (2?rA/ t T) cos( Afc) 

A 

E {Q 2i } w -yxDisinc (2tt Af { T) sin( A<fr ) (5) 

v2 

For an assumed orbital altitude of 1400 km and transmissions in the S-band, we expect a 
maximum phase acceleration due to Doppler of ~ 100. 2tt rad/s 2 . For T = 1 ms, we expect 
3tnc(27rA/i) > 0.97 so the factor can be safely dropped. The samples 104a,b are input 
to the microprocessor 56 which estimates the phase error and implements a loop filter to 
achieve desired phase- lock loop performance. Since the signal of equ. (1) has a common 
outer data sequence 86 on the in- phase and quadrature components, a simple Costas 
Loop discriminator approximates the phase error by multiplying the samples 104a, b. 

Sfa = I 2 iQ2i 
A 2 M 2 

E{Hi} = —j^-Afa (6) 
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From equ. (6), we see that the gain of the discriminator is K = A ^ • * n addition, the 
variance of the discriminator output can be computed 

°U = & {4.41} - & V*Q*) 

K 2 N 0 (. Nc 



A 2 T 



By monitoring the signal amplitudes of the correlator outputs 104a,b, and varying the 
loop control accordingly, a specific phase-lock loop transfer function, if, is maintained by 
the microprocessor 56. Generally, in selecting the bandwidth of the loop, a tradeoff must 
be made between rejecting thermal noise on the one hand, and tracking performance on 
the other. For loop transfer function H, the NCO phase error, which we cannot directly 
estimate by the technique described in section (VLB), has variance 



*4 



+ jj jT° S»(u) I 1 - H(jw) f <L> (8) 

where S n (0) is the power spectral density near the origin of the thermal noise. The 
expression assumes that the loop bandwidth, B L = ^ Jo°° I \ 2 du>, is much smaller 

than Be Since the correlator bandwidth, ^, is also much less than B c , we can safely 
approximate 5 n (0) = Tcr}^ From equ. (7), the phase error variance in equ. (8) becomes 



2 2N 0 B L f, . N 0 \ 



This thermal noise variance is not heavily dependent on there being the same data se- 
quence, {£>»}, on in-phase and quadrature signal components. For example, for a CDMA 
signal with QPSK data modulation, we would use a fourth-power loop, rather than the 
Costas Loop. It has been shown by Lindsey [16] that thermal noise variance could then 
be approximated: 1 + 3.8 (jfe) + 7.6 (^) 2 + 8 (jfe)*] . 

A phase-counter 57a-c keeps track of the absolute phase of the NCO 51a-c locked to each 
signal. The phase measurement of the phase counter 57a-c contains both a fractional and 
an integer component, where each integer refers to a full 2it phase cycle. There are multiple 
different means by which the microprocessor 56 can read the phase of the phase counters 
57a-c. Whichever method is employed, for each reading epoch, the time from reading 
the first phase counter tracking a signal to the last phase counter tracking a signal should 
span an interval less than a few psec in order to cause errors below the experimental noise 
floor. This specification is made concrete in section (VI. A). One method of satisfying this 
specification is illustrated in Fig. 6. The microprocessor 56 sequentially reads the phase of 
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each tracking module that is locked on a signal. This is achieved by means of a select signal 
58 which is input to the bus interface 59 of each of the phase counting assemblies 53a-d. 
Each select signal selects the output of one of the phase counters in one of the phase 
counting assemblies. The clock driving the phase 86 and select 58 buses should be fast 
enough to enable reading of all of the relevant phase counters in all of the phase counting 
assemblies in the specified interval. Another approach to the problem is illustrated in 
Fig. 7. In this configuration, a single latch signal is used to latch 64a-c the phase of each 
of the counters 62a-c in all of the phase- counting assemblies 53a-d simultaneously. The 
latched data 61a-c can then be read sequentially over a longer interval, by sending select 
signals 58 to the bus interface 63 for each of the relevant latches. Phase reading epochs, 
in which the phase data for all of the signals being tracked is read, occur at frequencies 
roughly between 1 and 100 Hz, depending on the requirements of the specific application. 

The frequency synthesizer subassembly 65 for the receiver 3,4 is driven by a single crystal 
oscillator 66. Many different techniques are known in the art for generating a series of 
frequencies from the crystal oscillator output 67. One method of generating the desired 
frequencies of a receiver path for the Globalstar LEO constellation is discussed in section 
(VIII). A real-time clock 68 is maintained continuously, and generates time information 
for the microprocessor 56 to aid in the initial acquisition of satellites. 

The microprocessor functions and memory for the user 4 and reference 3 receivers can be 
decomposed as illustrated in the block diagrams of Fig's 9 and 10. Each microprocessor 
consists of a CPU 73a,b, memory 70a,b,71a,b and a bus interface 74a,b through which 
to communicate with the rest of the receiver. The memory of the microprocessor stores 
the software routines 70 a, b. as well as the data 71 a, b necessary for the implementation 
of those routines. We will first consider the software routines 70a for the user receiver. 

70a. 1 - The initial acquisition algorithm performs a signal search routine to initially estab- 
lish phase-lock on satellite signals. This involves using the satellite ephemeris data 71a.l, 
the GPS signal structure data 71a.2, the LEO signal structure data 71a.3, possibly the 
location of ground uplink stations 71a.4, together with data 75a-d from the tracking 
assemblies 44a-d to control the tracking modules such that phase lock is established on 
the available satellite signals. Control commands are applied via the select /control bus 
76. The initial acquisition can also involve implementing frequency-locked loops as well 
as delay-locked loops, and is well- understood in the art (see, e.g. Dierendonck [17]). 

70a. 2 - The maintenance of phase-locked loops involves controlling the components of the 
tracking assemblies 44a-d so that phase lock is maintained on all relevant signals. The 
microprocessor implements a control law for closing the phase-locked loops, and possibly, 
delay- locked loops for each of the tracking modules in use. These control techniques are 
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well- understood in the art. 

70a.3 - Interpreting and demodulating data from the satellite downlinks involves using 
the data about the LEO 71a.3 and GPS signal structures 71a.2 to read and interpret 
the outputs of the tracking modules 75a-d. 

70a. 4 - The carrier phase measurement routine is the process whereby carrier phase data 
is read from the phase counters of each of the phase-counting assemblies 53a-d and 
interpreted to produce phase measurements that can be input into the data-reduction 
routines 70a. 6. 

70a.5 - The code phase measuring and positioning routine is the process whereby code 
phase data is read from the GPS tracking assembly 44a, corrected using the differential 
corrections received from the reference and processed to obtain a meter- level position 
estimate and clock offset estimate. 

70a. 6 - The position calculation routine involves the correction of carrier phase measure- 
ments as well as an estimation algorithm. The correction of carrier- phase measurements 
is the process whereby carrier-phase data is corrected for deterministic disturbances based 
on information received from the reference, satellite ephemeris information 71a. 1, possi- 
bly the location of ground uplink stations 71a. 4 and possibly error prediction data 71a.5. 
The nature of these corrections are discussed in sections (VI.A,VII.A,VII.B). The esti- 
mation algorithm is the key process by which data is processed to identify the integer 
cycle ambiguities for the GPS satellites as well as related parameters for the LEO satel- 
lites, and to subsequently position the user with centimeter-level precision relative to the 
reference receiver 3. The algorithm employs the satellite ephemeris information 71a.l, 
the GPS signal structure data 71a.2 the LEO signal structure data 71a,3, and possibly 
the location of ground uplink stations71a.4. The algorithm is discussed in detail in the 
section VLB. 

70a. 7 - Receiver autonomous integrity monitoring may be used by the receiver to indepen- 
dently check the validity of the position solution using the satellite ephemeris information. 
This technique is well-understood in the art and is described in more detail in section 
(VIILB). 

70a.8 - A phase velocity measurement routine may be employed to enhance the position 
estimation methodology. 

70a.9 - Filtering position estimates involves applying to the position data a digital filter, 
which takes into account known aspects of the users' motion, such as bandwidth con- 
straints, in order to generate more accurate, noise-free, position estimates. This could 
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also involve the use of Kalman-filtering techniques to combine the carrier-phase position 
estimates with data obtained from other sensors such as accelerometers and gyros. 

Most of the routines discussed also make use of assorted miscellaneous variables 71a.6 
We turn our attention now to the reference receiver, to discuss those routines which are 
not necessarily applicable for the user receiver. 

706.5 - The LEO clock calibration routine is used to identify frequency offsets of the LEO 
satellite oscillators. The algorithm makes use of the satellite ephemeris information 71b, 1, 
the GPS signal structure data 71b. 2, the LEO signal structure data 71b. 3 and possibly 
the location of ground uplink transmitters 71b. 4 to identify the frequency offset of the 
LEO downlink due to the long-term instability of the satellite oscillator. The technique 
is detailed in section (VII. C). 

706.6 - The code phase measurements and differential correction calculation involve read- 
ing the code phase measurements from the GPS tracking assembly 44a and combining 
this information with the known position data 71b. 6 for the reference antenna 17d, to 
calculate differential corrections for the user. 

706.7 - Transmitting data to the user is the process of coding and transmitting the data 
destined for the user receiver 4 in terms of the data communications protocol of the LRU 
8 that is employed. 

We have described the microprocessor operation assuming that all of the code is imple- 
mented on a single microprocessor. Another implementation might have multiple micro- 
processors in the receiver, each with specialized tasks. For example, a microprocessor for 
each tracking assembly might maintain signal tracking loops, while a separate micropro- 
cessor would be dedicated to the computation-intensive data reduction routines. Multiple 
permutations on this theme are possible. 

VL Precise Position Calculation 

VI. A. Details of the Differential Carrier Phase Measurement 

The scenarios in which centimeter- level positioning accuracy is required within a few 
minutes of system initialization include surveying, construction, precise control of land 
vehicles, as well as high-integrity tasks such as attitude determination and automatic 
landing of aircraft. The data reduction technique is similar in all cases, with simplifications 
for the context of a user receiver that is stationary in Earth- Based- Earth- Fixed (EBEF) 
coordinates. We will focus our attention on the more general case of a mobile user receiver 
4, which is driven by an oscillator 66 that is distinct from that of the reference receiver 
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3. The purpose of this section is to characterize the software and hardware upon which 
centimeter-level navigation accuracy depends. 

We denote the nominal carrier frequency of the satellite downlink as u a . The phase of the 
satellite frequency synthesizer at time t s can be described: 

rU 

^s{t s ) = u s t 8 + J Av s (a)da (10) 

where Au> 5 (£ 5 ) models the deviation from nominal frequency due to drift of the crystal 
oscillator onboard the satellite. We model this drift as a time- varying offset in the satellite 
clock, r s (t s ) = ~ J** Au>(a)da. Consequently, 

= u> 9 [t. + T,(t 9 )) (11) 

Similar clock offsets occur at the user and reference receivers. At true times t u and t r the 
user and reference receivers respectively record times t u + r v (t u ) and t r + T r {t r ). Therefore, 
the phase output from the crystal oscillator 66 of the user's receiver 4 at time t u is 

Vx(tu)="x[tu + T u {t u )] (12) 

where u>x is the nominal oscillator frequency. The frequency synthesizer 65 of the receiver 
generates the rf and if mixing signals by multiplying ^x(tu) by factors a r j and a,/ 
respectively. We denote the phase of the satellite signal emerging from the user receiver's 
LNA 22 at time t u as S u(K)- The output of the first mixer 26a-d, after bandpass filtering 
25a-d. has phase ^\(t u ) = *su(*u) — otrjwx [tu + T tt (t tt )]. We assume for simplicity of 
explanation that the scheme of Fig. 4a is employed. Hence, the second mixer 40, after 
filtering 41, generates a signal with phase ^(tu) = ^su(<u) — (<*»■/ + <^is) w x [tu + r u (t u )]. 
This phase is tracked by the phase- locked loop of the tracking module for satellite s, and is 
read from the corresponding phase counter by the microprocessor 56. Since the nominal 
satellite frequency is u; fi , the nominal offset frequency of the signal tracked by the PLL is 
u> 0 = u) B — (a r f + otif)u>x- The microprocessor 56 differences the phase it reads, ^P 2 (*u), 
with the phase component caused by this offset frequency, u>o- Since the microprocessor's 
measure of time interval is directly affected by crystal 66 instability, it calculates this 
phase component as u>o [t u + r u (t u )]. The resulting phase measurement is 

<MM = ™o [t u + r u (t u )\ - tf 2 (*«) - N 9U 2tt 

= u> s [t tt + r u (t u )) - V su (t u ) - N w 2r (13) 

where we have included an integer cycle phase ambiguity, 7V 5U , since the microprocessor's 
initial phase measurement is modulus 2tc. Consider now the phase of the incident satellite 
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signal, $ su . The signal is affected by phase disturbances on the satellite-to-user path, as 
well as frequency-dependent phase lags in the receiver. In addition, the phase measured 
depends on the position of the satellite, r a , at the time of transmission, rather than at 
the time of reception. Applying these factors to equ. (11), the phase of the signal from 
satellite s emerging from the LNA 22 of the user receiver 4 at time t u is: 



Pbu I tu I + VsT 9 Uti ~ j - Psu - Tlsu 

(14) 



where /x su is the frequency- dependent phase delay of the receiver and n au is the error due 
to ionospheric and tropospheric delay as well as thermal noise and imperfect carrier-phase 
tracking in the receiver. n su Is actually time- varying, but we ignore this time dependence 
for now. In order to represent the signal path length, we have denoted by p av {t 0 ) the 
distance from user's current position at time t u to the satellite position at transmission 
time to, p S u{to) = 1 r u (t„) - r a (t 0 ) I- Since we cannot know exactly the location from 
where a satellite transmitted if we do not know t 0 , the precise calculation of t 0 requires 
an infinite regression. However, one can make the simplifying approximation 



t„ - -P~ (tu - Ptu(tv "A « <« - -PsM 
C \ C } C 



(15) 



which generates a worst-case ranging error < 2 mm for satellites at 1400 km (the Global- 
star nominal orbital altitude). Therefore, from equ.'s (13) and (14) we estimate the user's 
phase measurement at time t u 



4> 9 u(tu) = ~P°* I " c J +Ws r u(K) - r s ^ u j 



-N su 27r+fx su +n su (16) 

The receiver assigns a timetag t to the measurement made at time t u . Since the user 
receiver does not know true time i u , the timetag will be effected by the clock offset of 
the receiver. The measurement must therefore be recast in terms of the clock of the 
user's receiver. If a user receiver makes code-phase measurements on the GPS signals, 
the clock offset in the receiver, T u (t u ), can be estimated to within 1/xsec. Once r u (t u ) 
has been estimated, two algorithmic approaches are possible. Firstly, one can use this 
estimate to select the times at which phase data is read from the phase counters in the 
receiver, so as to continually correct for the receiver's clock offset. This clock steering 
technique limits the magnitude of the timetag error, \t — t u \, to roughly Ifisec. Secondlj', 
one can use the estimate of T u (t u ) to actively correct for any errors which would arise in 
the differential phase measurement. Since these approaches are conceptually very similar, 
we will describe the latter approach in detail. Once the issues are identified and resolved, 
it will be clear to one skilled in the art how the algorithm varies for the former approach. 
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Assume the user receiver 4 estimates its clock offset to be f u using code-phase measure- 
ments. We assume for now that f u is time-independent since it need not be continually 
updated. We define Ar u (t u ) = r u (t u ) — t u . For phase data read at true time the user 
receiver assigns a timetag t = t u + Ar u (t u ) + T u , where T u is a residual error in sampling 
time, distinct from the clock bias, resulting from the digital logic's imperfect precision in 
reading phase at the particular instant of time identified by the microprocessor. Recasting 
the measurement in terms of timetag t: 

4>su(t u ) = ^Psu (t - Ar u (t u ) - T u - [t - Ar u (t u ) - T u )j 
+ UsT u (t - Ar u (t u ) - T u ) - N^Ztt + /isu 4- n au 
- w s t 3 (t - Ar u {t u ) -T u - ^p su (t - Ar u (t u ) - T u )) (17) 

Since \Ar n (t u ) + T u \ is of the order of a few /zsec, we may Taylor-expand to first order in 
Ar u (t v ) + T u , with negligible error in dropping the higher-order terms, 

4>su(t u ) « (18) 

tp- (* - + - ^^[Ar^o + rj 



Note that for highly unstable oscillators, one might include in this expansion the terms 
-uj s ^^[Ar u {t) + T u ] and ~Lo s ^^-[Ar u (t) + T u ], which could be incorporated into the 
estimation algorithm. However, these terms are negligible in most implementations of the 
invention. 

A completely similar approach to the reference receiver's phase measurement yields the 
expansion 

<M*0 w (19) 

^Psr (t - + »sT r (t) ~ ^^[Ar r (t) + T r ] 

-u, Ts (t - 2^£1) - N^tt + n sr 

The reference receiver's phase measurement is timetagged with a time t and transmitted 
(or otherwise communicated) to the user. The user matches the timetags on the data and 
performs a single difference, which we now index with the timetag t rather than a true 
time. The differential measurement is then 

4>*{t) = -M*r) (20) 

- *~ (' - "f) - "t" (' - ^ ) 
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+ [r, (t - ^ ) - r, (t - ^)] + u, 5 - r r (*)] 

^ dp»u(*) u 3 dp ar (t) f . , 

" --^l^(t) + rj + T -^[Ar r (t) + r T ] 

- (JV SU - A^ sr ) 2tt + (/x 5U - M*r) 4- (n stx - n 5r ) 



We may rewrite the terms involving and as 



[2 , _ rJ + /a^t) _ ^ [Aru(0 + rJ (21) 



or equivalently as 



+ ^a^o |7 , _ r<) + ^) _ ap^tYj ± {Ar (t) + j, j (22) 

We will adopt the representation shown in equ. (21). It will be clear to one skilled in the 
art how the issues we describe can be transferred to equ. (22). The 1 5< term in equ. (21) 
can lead to large errors and is incorporated into the estimation strategy. The 2 nd term in 
equ. (21) can be directly calculated, and subtracted from the differential measurement. 
If the user and reference receivers are implemented with similar digital logic, the term 
|T r — T u \ can be made less than .lfisec. Hence the 3 rd term can be ignored with distance- 
equivalent errors < lmrn. Now consider the 4 th term. For a satellite at 1400 Arm, and a 
stationary user receiver lOfcm from the reference, the term ( dp %^ — ^at^) < SOm/s. 
To ensure that the 4 t/l term produces a worst-case distance-equivalent error < 1mm we 
must have | Ar u (t u ) + T u | < 20/xsec, so that the term can be safely ignored. Consider the 
phase-reading scheme displayed in Fig. 6. The time interval required to read all of the 
active phase-counters of all of the active phase counter assemblies 53a-d contributes to 
the magnitude of the differentially uncompensated term lAr^t^-hT^I. Hence, for a static 
user, this time interval should be roughly < 18/xsec. However, for a mobile user receiver 
moving at, say, 250m/s, the term (^f^ — dp *Q t ^ ) < 300m/s, and the time period should 
be roughly < 3^,sec. The rate at which the estimate f u needs to be updated so that Ar w (t) 
remains small depends on the stability of the receiver oscillator 66. For example, for a 
long- term oscillator stability of 1 : 10 7 , updates every 2 minutes and every 20 seconds are 
sufficient for the static user and mobile user respectively. 
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The terms in equ. (20) relating to the satellite clock offset can be expanded to first order 



T . r, (t - « {pA t) - P „(0). (23) 



For a satellite oscillator with long-term frequency stability of 1 : 10 6 we expect to get 
as large as 10"" 6 . This could cause distance-equivalent errors as large as 1cm. In section 
(VII. C), we describe a technique for calibrating the frequency offset of the LEO satellite 
oscillator so that the first-order expansion of equ. (23) can be calculated and subtracted 
out of the phase measurement. 

Eliminating all terms from measurement equ. (20) which are either negligible, or actively 
subtracted out of the measurement, the estimate of the resulting measurement is 

«.co « >. i'-'^f)- > (< - + - (• - Mt) - *- m 

- (N, v - N sr ) 2tt + {il su - fi ar ) + (n su - n 8r ) (24) 

To clarify the estimation strategy, we redefine the components of this measurement as 
follows: 

T(i) = T.(i)-T r (t) (25) 

N, = N„-N„ (26) 
M. = P>u ~ P*t (27) 

where X s is the nominal wavelength of the satellite carrier. Multiplying equation (24) by 
|j to convert phase to a distance measurement, we then have 



».w = p«u (* - a^i) - P „ (* - a^i) . (29) 

+ (l - i 2 "^) cr(t) - N s \, + t»> + n » 
VLB. Estimation Strategy 

If one attempts to estimate all of the integer components, {iV*}, as well as user position, 
r„(t) and clock biases r(t), the complete set of parameters would be almost unobservable. 
The resulting estimation matrix would be poorly conditioned and highly susceptible to 
measurement noise n a . Consequently, we select one of the Navstar satellites, say satellite 
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1, to be a reference satellite for differencing. We make an initial approximation of the 
associated integer, N\ using equ. (24), based on our estimate of position and clock offset 
using code-phase measurements. Then, we adjust the measurement 

Vs(t) = y.(t) + Ni (30) 
and we redefine the parameters we seek to estimate as follows: 

*• - - £) - > - & 

?(<) - <•(<) - 7- (th - M " §) (3D 

For the Navstar satellites, the new parameter N s simply reduces to N a — since all GPS 
satellites are transmitting similar frequencies. For the LEO satellites, the phase-delays 
do not cancel and the parameters N 3 cannot be regarded as having integer values. If 
| iVi — A^t — ^ | <: 200 cycles, the redefinition (31) leads to a maximum distance-equivalent 
error of < 1mm for satellites at 1400fcm. The measurement can then be approximated: 



!/.(0 - Psu 



+ U - cf(t) - N t \s + n t (32) 

where Ni is by definition 0. 

The set of time-dependent parameters which we seek to estimate is 

0(0= [r„(t) T cf(t)] T (33) 



We create an observation matrix for these parameters based on our estimates of the 
line-of-sight vectors to the satellites 

P.u (t - 1 dp n (t) ' 



Mi) = 



Pm (t-*-®) c at 



(34) 



For a current set of estimates, 0(0 and 7V a , we can construct an estimate of our prediction 
error for satellite s: 

Ay.(t) = p. u (t-5^)- P „( t -rf!) 



V c at ) 



+ 1 ~- ri ^r / )cf(t)-N 3 \,-y.{t) (35) 
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Estimation matrices and prediction errors for all visible satellites are stacked into com- 
bined matrices, 



H(t) = 









h 2 (t) 

• 
• 
* 


AY(t) = 


Ay 2 (*) 

• 
• 
• 


h s (t) 




. Ay s (0 . 



(36) 



where 5 is the total number of satellites visible. The batch measurement equation, involv- 
ing measurements from U through ijv and the batch parameter update, A0, can then be 
expressed as 

AY = HA0 + V (37) 



with the matrix structures 



AY = 



AY(M 
AY(< 2 ) 

AY(* W ) 



H = 



A© = 



H(ti) 0 
0 H(* 2 ) 0 



A0(tO 
A©(t 2 ) 

■ 

A0(tjv) 
AiV 2 



AN S 
—A 
— A 



H(* N ) -A 



(38) 



where 



A = 



A 2 0 
0 A 3 
0 
0 



0 
0 
0 
0 

As 



(39) 



The disturbance matrix V contains errors due to each satellite's measurement noise n a (t) 
- which we may reasonably assume is uncorrelated with distribution n s (t) ~ N(0^a^ 8 ) - as 
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well as ephemeris errors e s (t) due to imperfect knowledge of the satellite's position which 
affects calculation of p su and p sr . Combining these two noise sources, the disturbance 
matrix has the form: 

n x (U) + e^ti) 



V = 



n s (U) + e s (t 2 ) 
ni(t 2 ) + ei(t 2 ) 

• 
• 

ns(t 2 ) + e s (t 2 ) 



(40) 



It should be noted that the matrix structures and parameters described above can be 
altered if the user receiver 4 is static relative to the reference receiver 3. This situation 
pertains, for example, in surveying applications, or any problem where a vehicle can 
remain stationary until good integer estimates are available. In such scenarios we need 
only estimate the 3 coordinates of r u (*i). Given an estimate, r u (ti) of the initial position, 
our estimate of the position at time t n is simply G(t n — <i)f tl (ii), where G is a rotation 
about the z-axis in EBEF coordinates, which accounts for the earth's rotation between 
time < n and t\. To account for the reduction in the number of parameters, we define for 
each time t two separate stacked observation matrices: 



Hi(t) = 




H,(t) = 



1 - 



1 - 



c dt 



I dvsu(t) 
c dt 



(41) 



We then structure the batch estimation matrix, H, as 



H = 



Hi(*i) H 2 (h) 0 ••• 

H,{« a )G(* 3 -t,) 0 H 2 (tj) 0 

; o 

Hi(t w )G(t w -ti) 0 



—A 
—A 



H2(<n) —A 



(42) 



and the batch parameter update matrix, A© as 



S» ii 



A© = [Ar x (*0 Af y (*i) Af z (ti) cAf(iO ... cAf(t N ) AW 2 ... AN S ] 



(43) 



26 



p 



♦ * I 

WO 98/43372 



PCT/US98/06042 



and we may proceed with a batch measurement equation as in (37) above. 

Once the estimation problem has been formulated in the manner of equ. (37) the solution 
can be well-conditioned. The good conditioning is due to the geometric diversity resulting 
from the motion of the LEO satellites. This geometric diversity decreases the condition 
number of the estimation matrix H 



k{h) = ^ox(Hl (44 ) 



mtn 



where cr mtn (H) and c7 max (H) are the minimum and maximum singular values of H. The 
condition number indicates the sensitivity of the parameter solution to disturbances V as 
well as to errors in the estimation matrix, 6H. This concept can be made more concrete 
by considering the || H2 norm of the parameter estimate errors for a simple least-squares 
parameter solution. Imagine A0* is the parameter update solution of the least-squares 
problem with these error sources removed: 



while A0£,s is the actual least-squares solution found 



A0 LS = arg min ||HA© - AY|| 2 (46) 



It can be shown (see Golub [IS]) that if 



e = max 



Iw-iiAYjur^ (47) 



and 



JIHAe^AYIb 
P ||AY|| 2 K ™> 



then 



\\A@ LS - A0.|| 2 

\\*®LS\\2 



5 £K(H) (^8) + tan </? )*(H)) ( 49 ) 



For one skilled in the art, equ. (49) indicates how small /c(H) should be for a given t in 
order to achieve a desired level of accuracy in the parameter estimates. 
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VI. C. Mathematical methods for solving the estimation problem 

Many different mathematical methods can yield a solution to the problem posed in equ. 
(37). A method can be chosen depending on the processing power available in the receiver 
and the requirements of the specific application. Some of these different approaches are 
discussed below. These in no way represent a complete set of techniques; rather they 
highlight some key methods. The parameter solution found by any of these methods 
would be well-conditioned due to the geometric diversity achieved by the LEOS, based 
on the reasoning of equ. (49). 

VI.C.l Maximum Likelihood Update 

We assume now that our estimated parameter set 0 is near enough the true parameter 
solution that there are negligible errors due to higher order terms, incurred in lineariz- 
ing the measurement equation to derive (37). Based thereupon, we seek the maximum 
likelihood update: 

A® ml = axg max Prob { AY | AG)) (50) 

This ML update requires knowledge of the measurement covariance matrix. Since the 
ephemeris errors are strongly correlated, the covariance matrix C = .E jw r j has non- 
diagonal structure. Over the course of 5 minutes of tracking, one can roughly model the 
satellite position generated from the ephemeris data as 

r.(*) = r,(t) + Ar. (51) 

where Ar s is a constant offset, modeling the difference between a satellite's true position 
and the position estimate based on the ephemeris data. We describe this Offset vector in 
terms of normally distributed components: 

Ar=[Ax Ay Az] T , Ax,Ay,Az - w(o,<r*). (52) 



This error in the ephemeris data would result in an ephemeris disturbance - see equ. (40): 

C.(*) = (Psu(t) - Psr(t)) ~ (Psu(t) ~ Psr(t)) (53) 

ArsjtfrUt) Ar 8 (t) T r sr (t) 

r S u(t) r sr {t) K } 

where the approximation is achieved with a first-order expansion, assuming that || Ar s || 2 <C 
r su, r sr- Using this first- order approximation, we find the second moment of the ephemeris 
disturbance statistics: 
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£{e 5 (ti)e 5 (< 2 )} = a\ 



r n (ti) T r tu (t 2 ) r„{U) T r tT {t 2 ) r„(ti) r r„(t 2 ) r„(U) T r su (t 2 ) 
r s „(ii)r sti (t 2 ) r t , T (t 1 )r„(t 2 ) r„(ti)r„(t 2 ) r„(ti)r s »(t 2 ) 



(55) 



Consequently, the batch covariance matrix C will have structure: 



C = 



C n + C e (z,) C.fr.ti) C.(t u t z ) 
C.(*2,*i) C n + C e (i 2 ) C e (i 2 ,* 3 ) 
C«(t 3 ,*i) C e (t 3 ,i 2 ) C„ + C e (t 3 ) 



(56) 



where 



C.(t<,<j) = 



0 



0 



(57) 



C n = 



Given matrix C, the ML parameter update is: 



«2s 



A©ml = (^C^H) 1 H r C _1 AY 



(58) 



(59) 



In essence, this iterative estimation strategy is the Gauss-Newton technique, where we 
have pre- multiplied the batch estimation equation (37) by the whitening matrix C"2 to 
achieve the ML update. Since the measurements in equ. (32) are only mildly nonlinear 
in r„, 2 — 3 iterations are sufficient to converge to the experimental noise floor. 

VI. C. 2 Least-Squares Batch Solution via Choletsky Factorization 

The technique described above for solving the ML estimation problem requires 0(N 3 ) 
flops. Many techniques exist for reducing the computation time required by exploiting 
the sparse structure of the measurement matrix H. One such technique is discussed in 
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this section. In order to preserve sparseness. we ignore the off-diagonal terms of the batch 
covariance matrix of equ. (56), to obtain a diagonal matrix C. We may pre-multiply the 

batch measui'ement equation by the diagonal scaling matrix (c) 2 and then solve the 

least-squares problem. Since (c) 2 does not change the block structure of H, we will not 

explicitly show the premultiplication by (c) 2 , or equivalently assume that C is simply 
the identity, Ins- We may then solve the least-squares problem of equ. (46) by solving 



H T HA© = H T AY 



(60) 



The matrix A = H T H has the block structure 



Am 
0 



0 

A 2 ,2 

0 



0 

A 3i3 



Atv+i.i Ayv+i^ Ayv+1,3 



Ai,N+l 
A 2 ,AM-l 

A 3> N+1 

« 
• 
• 

Ajv+i^+j 



(61) 



where the submatrices {A t| ;}, i < N + 1 are 4x4, { Ayy+i,,} are (S — 1) x 4 and A^/ + i t ^+i 
is (S — 1) x (5 — 1). We seek a lower triangular matrix, L, such that LL r = A, with 
structure 



L = 



0 



0 

^2,2 

0 



0 



0 



(62) 



where the submatrices {£,,<}, i < + 1 are 4 x4 lower triangular, {Lw+ 1|t } are (S-l)x4 
and L/v+i,n+i is (5—1) x (5— 1) lower triangular. This matrix can be found via Choletsky 
Block Factorization (see Golub [18]), which is achieved with the following algorithm: 



for j = l,7V + l 
for i = j, N + 1 

S = Aij - EUi L^LT, 
if i = j 

compute by Choletsky factorization L JtJ s.t. Ltjjhjj 



= S 
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else 

solve LijLjj 
end 
end 



= S for L tJ 



end 



Once L is found, A0 can be found by block back-substitution. 
VI.C.3 Iterative Information Smoothing 



Despite the computational efficiency, sparse matrix batch algorithms are difficult to ad- 
minister when satellites are coming in/out of view, or cycle slips occur, while data is 
being stacked for processing. Information smoothing might be selected for its flexibility 
in dealing with such situations. In essence, the information smoother passes the linear 
Kalman filter forward and backwards over the data before updating the parameter set. 
We describe the set of parameters as 



AX(t) = 



A0(t) 
AN 



(63) 



where AN = [AJV 2 . . . AJvs] T . We model the evolution of parameters between phase 
measurements as a Gauss-Markov process 



AX(/ t ) = AX^-O+W^) 



(64) 



where 



E{w{t k )W(t k f} = Q(t k ) 



= lim 



C tv — *00 



l40u, 2 0 



2 

u 

0 0 



(65) 



This process is based upon the idea that the integers are constant for all phase measure- 
ments and that no assumption is made about the user's motion between measurements. 
The phase prediction error of equ. (36) we model with a linear approximation: 



AY(t k ) = H(t k )AX{t k ) + V(r fc ) 



(66) 
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where 



E {V{t k )V(tk) T } = R(h) = 



0 



0 



(67) 



Denoting the covariance of our estimate AX(ffc) as P(**)> the Kalman filtering equations 
for this system are (see Kalman [19]): 

AX(t 4 ) = AX(ijk-i) + P(tjb)H(tfc) T R(^)~ 1 ^AY(^) — HffJAXft^i)) (68) 



If phase-lock is attained on a new satellite 5 during parameter estimation, the initial 
covariance of the parameter estimate AN a is very large; this can cause computational 
difficulties. Consequently, we use instead an information form of the Kalman filter (see 
Pervan [10] and Stengel [20] where we define an information matrix 



8(0 = P(0 



S&(t) S e #(t) 



(69) 



and an information vector Z(t) = S(/)AX(t). Applying these definitions to equations 
(65) and (68), it is straightforward to show that the update equations become: 



S(t t ) = 



0 0 

0 S N (*fc-i) - S QW (r k _ 1 ) T S e (i Jt _ 1 )- 1 S e //(<it-i) 



+ H(i fc ) J R(**)~ 1 H(* fc ) 



Z(t k ) = 



_ Z N {t k -i) -S e /K**-i) T Se(*k-i) _1 Ze(ijfe-i) 



+ H( lfc ) i R(<0" 1 AY(i fc ) 



(70) 



In order to emulate the batch solution, the filter is passed forward and backward over the 
data. For the backward pass, we simply interchange the tk~\ and tk in equ. (70), and 
start with initial conditions S(t^) = 0 and Z(*at) = 0. The S(t*) and Z(tk) from forward 
and backward passes are linearly combined 



S(i fc ) = S(t fc ) F + s(t fc ) B 

Z[t k ) = Z(t k f + Z{t k ) B 



(71) 
(72) 
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and the parameter updates are then found according to AX(tfc) = S(<jt)~ Z(tfc). The 
integer updates are found from the relevant vectors in AX(*i) or AX(iat). This process 
is repeated until the elements of AX(tfc) are negligibly small. 

VI. C. 4 Solving the Max. Likelihood Problem with Integer Parameters 

All of the techniques discussed thus far treat the integer parameters as real numbers in the 

estimation strategy. In this section, we discuss a technique for solving the least squares 

problem assuming an integer parameter set (see Hassibi [21]). Assume that our current 

parameter estimates are close to the ML solution, so that no iteration is necessary and 

we drop the A notation. Consider the batch measurement of equ.(37) rewritten as 

• ■ 

Y = Hfl0 + H z z 4- V (73) 

where He and H z are the estimation matrices for the real and integer parameters respec- 
tively. 9 is the matrix of real parameters, z is the q xl matrix of integer parameters, and 
V is the disturbance matrix with covariance C. Imagine equ. (73) is premultiplied by the 
whitening matrix, C~2 ? and the least-squares solution is found, with the approximation 
that z contains real elements. The resulting estimate of z can be shown to be normally 
distributed, z ~ 7V(z, *) where 

* = [hJc-'H, - HfC- 1 ^ (H^C- 1 ^)' 1 HjC-'U z y l (74) 

We can consider the least-squares solution z to be generated by z = z + u where u is the 
disturbance term. Multiplying by the whitening matrix G = , we have z = Gz + u 
where u ~ iV(0, I 9 ). The maximum likelihood estimate of z is then 

Zml = arg min ||z - Gz|| 2 (75) 
xez« xl 



G forms the basis, or generation matrix, of a lattice, L(G) = {Gz | z G Z 9 * 1 } which 
is conceptually illustrated in Fig. 11 for q = 2. Given the lattice L(G), 108, we can 
find an integer- valued matrix P which maps F : Z 9 * 1 -+ Z 9 * 1 such that \det(F)\ = 1 
and L(GF) = L(G), 108. Given the matrix G, the algorithm of Lenstra, Lenstra and 
Lovasz (LLL) (see Hassibi [21]) can be used to find an F such that the generation matrix 
G = GF has two key properties: 

•The columns of G, {g,}, are almost orthogonal, namely gy = £{=i Mjig'i Pii = *i < 
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\ for i •£ j where {g t *} are the columns of G*, the Gramm- Schmidt orthogonalization of 
G. 

•The || 1| 2 norms of the collumns of G are bounded, namely ||gi|| 2 * ■ ■ llg^lb < 2*<*- 1 >/ 4 |det(G)| 
and \\gih < 2^-Wrf\det(G)\. 

Before actually determining the integer solution, we desire a lower bound on the Prob{zML = 
z). To determine this, we seek the radius 107, of the largest ball 106 that fits with- 
in a Voronoi cell 109 of lattice L 108. If G is orthogonaiized using the Gram-Schmidt 
algorithm to G* = [gjgj . - . gj] , a lower bound on <* mtn can be found: 

d mtn >min{||gl || 2 ,|| g ; |! 2 ... || g ;|| 3 } (76) 



By redefining the generation matrix G GF, we tighten the bound of equ. (76). The 
sum-of-squares of q independent normally distributed unit variance random variables is 
a x 2 distribution with q degrees of freedom. We denote by F x i{x 2 \q) the cumulative 
distribution function of a \' 2 random variable with q degrees of freedom. Once we have a 
bound on c/ m in, d y we can find a lower bound on the probability of correct integer selection: 

Prob {zml = z} > F X 2 qj (77) 

We return now to the integer least-squares estimation problem of equ. (75). Replacing 
G by G, the expression can be rewritten as 

zml = arg min (z - z) T P~ l (z - z) (78) 

where z = G _1 z and P = (G T G) _1 . If P is diagonal, then the expression -becomes 

*ml = arg min £ ' p (79) 

in which case we can simply find the integers by rounding: zmu = F 2 «J • Since G is almost 
orthogonal, we can use the rounding of z as an initial suboptimal estimate of zmIo z *t*&- 
Since we have the lower bound ^y* from equ. (76), we know that if ||z - Gz, U 6|| 2 < 
then z, u 6 is the global minimum zml- Consequently, an efficient algorithm for attaining 
the global minimum is as follows (see Hassibi [21]): 

* Perform the LLL on the collumns of G to generate a unimodular matrix F and create 
a new lattice generator matrix G = GF which is almost orthogonal. 

* z 8ub *- F[G _1 zJ. If \\z - Gz^ fc || 2 < then z M l = z sub and stop. 
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* Using the algorithm of Babai [22], find a new z Jtl 6, such that 
z - Gz^ = ELi A t -g? with | Aj |< 1 

* If ||z - Gz 5U 6||2 < then z M L = z su6 and stop. 

* r 4- ]|z - Gz, u6 || 2 and z* *- F^z.^ 

* while 1 

* Search for an integral point inside an ellipsoid ||z — Gz|| 2 < r. 
If no such point exists Zml Fz* and stop. 

* Let z* be the integral point found in the previous step, r *— 1| z — Gz* || 2 . 

* If r < then z W l <— Fz* and stop. 

* end. 

Once the parameters {N*} are identified using any of the techniques discussed above, 
they can be regarded as constant biases. The receiver uses these estimates in equ. (35) 
to construct prediction error estimates with which to update position estimates with 
centimeter-level precision in real time. One straightforward method of doing this is 

A®{t k ) = (HCtO'R^'HftO)" 1 U(t k ) T R(t k y l AY{t k ) (80) 

VII. Implementing the Invention with Different Communication 

Configurations 

VILA. Satellites with Multiple- Beam Downlinks 

In section (VI), it was assumed that each satellite footprint was constructe'd with a single 
beam so that continuous carrier-phase data could be accumulated over the duration of 
tracking. However, rather than transmit a single beam for the satellite's footprint, as in 
the GPS case, many LEOS transmit multiple beams. As shown in Fig. 12 which roughly 
resembles the Globalstar downlink, each beam 79 a- d covers a portion of the satellite 
footprint 78. We may assume that each beam 79a-d is modulated differently, and has a 
phase offset relative to adjacent beams. Imagine a receiver moves from beam a to beam 
b 80. There arises an interval of a few seconds between reaching the 3cf5-down point 
of beam b 81 and the 3di?-down point of beam a 82. During this interval, the receiver 
simultaneously tracks both beams. Before handover from beam a to beam 6 at time to, 
the receiver calculates the phase difference between the two beams, 4> a {t 0 ) — <£&(*o). The 
receiver then adds this difference to the phase measured on beam b at some later time t\. 
The sum then becomes <f> a (to) + [&(ii) — $a(*o)]> which is corrected for any offset between 
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the beams. 

VII. B. Compensating for the effects of a bent-pipe communication payload 

We have assumed in section (VI) that the signal arises onboard the satellite. However, 
for a bent-pipe communication payload, the satellite actually downconverts and then re- 
transmits signals received from ground-based uplink stations. This does not change the 
conceptual approach, but simply requires that additional terms be taken into account in 
the estimation strategy. An example of a bent-pipe system configuration, such as that of 
Globalstar, is illustrates in Fig. 13. Consider the phase of the signal 84a, b generated at 
the ground terminal 83a,b: 

»l(«f) = »9 1*9 + T,{t,)] (81) 

We will ignore r 3 since its effect on the differential phase measurement is negligible. Using 
similar reasoning to that of section (VI. A), the phase of the uplink signal 84a,b incident 
at the satellite 2a,b at time t s can be described: 

<MO = (*. ~ \pM) + v 9 * (82) 

Where p gs (*s) 15 tne distance from the ground terminal to the satellite at time t s \ the 
expression t s — ^p g8 (ts) describes the time at which the transmission was made; and v g8 
represents all the phase disturbances on the path from ground to satellite, which will 
almost cancel out in the differential measurement. The satellite downconverter mixes the 
incident signal with another at frequency u> g — u; s . The phase of the down- converted signal 
which the satellite transmits is then 



*.(*.) = W. "[*. + Ts(U)) - — P 9 s{ts) ~ »gTs(ts) + *>,. (83) 

C 

The phase of the satellite signal output from the LNA 22 of the user's receiver 4 at time 
t u is then approximately 

(84) 

where v su contains all the phase disturbances on the satellite-to-user path. Consequently, 
we find the user's phase measurement corresponding to equ. (16) 

4>su(tu) = a; 5 r u (t u ) + — psu Uv ^ J + K - ^s)r 8 U« 1 
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+ ^P 9 . U. - Es^j - Ug , N„2* (85) 

Similarly, the phase measurement made at the reference at time < r is 

4,„ { t T ) = W .r r (M + f P .(t,-^)+K--.)r.(^-^) 

+ % ff5 ^-^^)-^-^-^V 5r 2 W (86) 

As discussed in section (VI. A), the user matches the tags on the measurements, and 
performs the single difference: 

M*) = Mtu) - MM (87) 

Using a Taylor expansion and discarding all insignificant high-order terms, as in section 
(VI. A), we can recast the measurement in terms of the assigned timetags. The resultant 
representation of the measurement, corresponding to equ. (20) is: 

c ot c . ot 

+ u,. [r u (0 - r,(*)] + K - «.) (' " ^) - r. (* - ^)] 

(88) 

We will consider those terms in the expression which have been introduced by the bent- 
pipe architecture and cause equ. (88) to differ from equ. (20). Consider the terms 



_ [ T . (« - ^2) - r. (< - *f>)] (89) 



involving the instabilities of the satellite oscillators. All that has changed from equ. (20) is 
the multiplication factor, which is now (u g —lo s ) instead of — u> a . Therefore, the treatment 
of this term is similar to that of equ. (23). 

Consider the terms in equ. (88) involving We rewrite these terms as: 



at 

^K(t)-r r (t)] + ^[f.-^] 

<fg dPgsjt) ^ _ T , _ U g dp gs {t) f p 8U {t) _ Psrjt) . ^ 

c dt c dt \ c c 1 



) 
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The 2 nd can be directly calculated if the position of the ground station and the satellite 
ephemerides are known; hence it can be actively subtracted out of the measurement. The 
3 r d term can be ignored, using similar arguments to those of section (VLA). The 4 th 
term is directly calculable if the receiver knows the satellite ephemeris, the location of the 
ground uplink, and the rough position of the user using code-phase measurements. Con- 
sequently, the fourth term is also calculated and actively subtracted out of the differential 
measurement. Only the I st term must be directly estimated. 

Eliminating all terms from equ. (88) which are either negligible, or actively subtracted 
out, the estimate of the resulting measurement corresponding to equ. (24) is 

+ [ Tv (t)-T r (t)} 

- (N sv - N„) 1-k + (/x ail - n„) + (n 5U - n sr ) (91) 

Following the steps outlined in section (VLA) and section (VLB), we find the measurement 
approximation corresponding to equ. (32): 

y*(t) = Psm It ^ 1 - p S T l< ^ 1 + 

(l - \^ - ^l ^) - NsK + n s (92) 



The observation matrix for the time- varying parameters, corresponding to equ. (34), 
becomes 

P. u (t - ^) 1 dj>»(t) co, d Ps ,{t) 

p„ (t - E^l) c dt u> s c dt 



h.(t) = 



(93) 



r 

Similarly, the estimate of our prediction error, corresponding to equ. (35) becomes 

A»(t) = P~ (* " ^) " Psr (t - (94) 

+ i 1 -7-9r"-^-ar-j cr(<) - 7V ' Ai - y ^ 

and we may proceed with the estimation as described in section (VLB). 

VII. C. Unstable Oscillators: Calibrating the LEO Oscillator using Navstar Satellites 

In this section, we describe how the frequency offset, or clock offset rate of the LEO 
oscillator can be calibrated using a GPS signal. This algorithm is only necessary for 
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oscillators that have long-term frequency stabilities of 1 : 10 6 or worse. Some of the 
mathematical steps are closely related to those described in detail above, so these stages 
have been left out of the explanation. The technique described here has been designed to 
be implemented completely with software, and requires no additional front-end hardware 
in the receiver. We will assume a bent-pipe communication architecture for generality; 
the additional terms can simply be dropped for simpler systems. 

We can describe the phase measurements made for a bent-pipe LEO satellite L, 2a-c, 
and a Navstar satellite N la-d at the reference receiver 3 as: 

4>Lr(tr) = [U + T r (t r )] - * L r(*r) - N^IT (95) 
<f>Nr{tr) = "AT [tr + T r {t r )] - *Nr(*r) - N Nr 2lt (96) 

In order to cancel out the error due to the receiver's oscillator 66 drift, the microprocessor 
56 performs a weighted difference between the phase of the two satellite signals to find a 
calibration phase 



4>c{tr) = <t>Lr(tr) ~ —<f>Nr{U) (97) 

The incident phase from each of the the satellites may be described: 



*Lr(*r) = 10 L 



(98) 



*Nr(*r) = U N lt T - -P/V r U r - J I + U)ffT N I t T - J - jlNr n Nr 

where the subscript g refers to the ground uplink station 83a,b, discussed in section 
(VII.B). The resultant expression for the weighted difference is found to be: 



JL /„ \ (4 PNr{tr) \ / PNr{tr) \ t*>L ( PLr(tr) \ 

MU) = «LTN^*r —)-T PNr \ tT 7~ )+T PLr \ tT 7~) 

+ -fP 9 L Ur " 1 + fog — UL)T L Ur " I 

- N^w + —N Nt 2k + mit -fiNr + n Lr n Nr (99) 
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The 2 nd , 3 rd and 4 t/l terms are directly calculated and subtracted from the measurement 
by the reference receiver which knows the location of the ground uplink station, as well as 
that of the LEO and Navstar satellite. After subtracting out these terms, the calibration 
measurement becomes: 



- N Lr 2* + — N Nr 2ir + fi Lr ti Nr + n Lr n Nr (100) 

For the purpose of this calibration, we consider the LEO's clock offset, ti,, to be a linear 
function of time as a result of a frequency offset in the satellite oscillator. Since the GPS 
clock is atomic, we regard as constant. The reference receiver calculates the change in 
calibration phase A<£ c over an interval of roughly one second, At, to find 

— ^ {Ug - UL )— (101) 

from which ^ can be calculated with good accuracy. 
VII. D. LEO Satellites using TDM A Downlinks 

It should be noted that the fundamental technique of augmenting GPS with LEOS for 
geometric diversity is equally applicable to TDMA downlinks, where the LEO satellite 
signals 6a,b arrive in bursts of a few psec. The time from the start of one burst to 
the start of the next is termed the scan period, T s . The time duration of each burst 
is termed the receive time, T r . It has been demonstrated by Cohen [5] that continuous 
carrier-phase tracking of GPS C/A code-type signals can be achieved with high integrity 
for T r = 2msec, and T 8 = 12msec. Although LEOS signals show more Doppler shift 
than GPS signals, it is well known in the art that a 3 rd order phase-locked loop can 
be implemented to maintain a running estimate of phase 0, phase rate -g£, and phase 
acceleration Hence, the change in phase due to satellite motion can be estimated 

over T 6 to maintain phase lock. Two fundamental limitations on the technique exist. The 
first concerns the stability of the satellite and receiver oscillators required to guarantee 
that cycle slips do not occur between bursts. If A is the Allen variance of the limiting 
oscillator, we require 

« 1 (102) 

27T 

For Ka band downlinks for which ~^ « 30 GHz, and T s ~ 25msec, we require i4<1.3x 
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10" 9 . The second limitation involves aliasing due to the dynamics of the receiver platform. 
For T s = 25msec, for example, the highest frequency component of the platform dynamics 
should not exceed half of the corresponding sampling frequency, or 20Hz. Neither of these 
constraints is restrictive for the majority of applications. 

VILE. Non~GPS Navigation Signals 

In our discussion of the carrier-phase positioning algorithm, we have assumed that LEOS 
are used to augment the Navstar GPS satellite fleet. It should be noted that Navstar 
GPS satellites have utility to the extent that they provide: 

•Additional carrier-phase signal sources, which in essence provide additional equations so 
that the parameter-estimation problem is over-determined. 

•Code-phase navigation signals which allow for correlation of the user and reference receiver 
clocks to better than ~ Iftsec and also enable the receivers to achieve initial position 
estimates accurate to the meter-level. 

Other navigation satellites exist, and are planned, which could also fulfill both of these 
functions. Such systems include Russia's Glonass, and Europe's proposed Satellite Civil- 
ian Navigation System which will include a LEO segment with highly stable downlink 
frequencies. Our system of exploiting LEO satellites for geometric diversity is equally 
applicable to these other navigation satellites. For the general case of a mobile user, 4 
such navigation signals should be available. If another means is used to initially synchro- 
nize the reference and receiver clocks, specialized code phase navigation signals become 
unnecessary. 

The technique described in section(VI) can resolve integer cycle ambiguities on carrier- 
phase signals for a mobile user so long as a total of 5 satellites are in view. Ideally, to 
rapidly constrain all the degrees of freedom in the positioning problem, 2 or more of these 
satellites should be in Low Earth Orbit. 

VII. F. Non- Differential Position Estimates 

The essential technique of augmenting GPS with LEOS for geometric diversity is equally 
applicable to the non- differential setting. To examine how a user receiver might proceed 
with non-differential position estimation, consider the first-order expansion of the phase 
measurement in equ. (18) 

MU) « ^(t-^) +Wjru ( t )-^ Mt)+ T„] 
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— u) 3 r s (t — P — N su 2ir + /x 5U + n dtl 



Since there is no need to initially estimate and compensate for clock offset, f u has been 
dropped from the expression. The term ^ a P^(0 r u can De ignored with distance-equivalent 
error of < 1mm for |T U | < 0.1 fisec and satellites at 1400 krn. Consider now r. 
For highly stable satellite oscillators, such as atomic oscillators, the user can be conveyed 
the necessary information to model the clock term. Alternatively, a highly stable clock 
can be modeled linearly as r 8 (t) = r so + r sl t. Using the clock-calibration technique 
described in section (VII. C), the user may directly estimate and subtract out the term 
(t — Ei ^)t s1 . Using Ionospheric/Tropospheric models (see Parkinson [23]) and/or dual 
frequency ionospheric calibration, the user receiver could estimate and largely subtract 
out the error terms n su , leaving a residual An su . Converting the resultant estimate to 
distance-equivalent form, 

VM = P,u (t» - + (l - - c ^^) «■-(*.) - - KN, U + ^ 5U + An 

(103) 

Using one of the Navstar satellites, say satellite 1, as a reference satellite for differencing, 
we make similar redefinitions to those described in section (VLB) 

y*u(tu) = y 8 (t u ) + Wi u 

f v (t u ) = T u (i u ) - ^ (N iu - Ny u - ^) (104) 

For the Navstar satellites, for which t 90 = 0, the new parameter N su remains integer- 
valued. The measurement may then be described 



VM = Psu (tu - J ^^j + (l - y 2 ^) cfn{t n ) - \,N SU + An tu 



(105) 



One skilled in the art could then proceed with a similar estimation strategy to that 
described in section (VLB), with the positioning accuracy depending primarily on the 
magnitude of An 5tt . 

VIII. Example Implementation of the Invention: Augmenting GPS with 

THE GLOBALSTAR TELECOMMUNICATIONS CONSTELLATION 
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VIII. A. The necessary criteria for cm-level positioning 

In order for a LEO constellation to rapidly resolve cycle ambiguities with integrity for a 
mobile user, the following criteria should be fulfilled: 

•There should ideally be 2 or more LEO satellites available for tracking. 
•A carrier signal should be traceable for a time period of a few minutes. 
•The satellite ephemerides should be known to good accuracy. 
•The SNR ratios should be sufficient for accurate carrier phase estimation. 

All of these criteria are fulfilled by the Globalstar Constellation. Carrier phase from one 
satellite can be tracked for several minutes at a time. In addition, the constellation has 
GPS sensors onboard so the ephemerides can be estimated to < 20m rms (see Yunck 
[24]) Fig. 14 describes the percentage availability of the satellites at different latitudes. 
Note that there are always 2 satellites available above 10 deg elevation over the continental 
United States. 

VII LB. Integrity with RAIM 

In all availability and performance analyses, we assume that the LEOS are functional, 
and that no cycle slips occur over the tracking duration. For high-integrity applications, 
this assumption cannot be made, and position solutions would be independently validated 
via receiver autonomous integrity monitoring (RAIM), (see Parkinson [25]). In essence, 
the RAIM algorithm checks if the residual of the least-squares solution at each time 
t, || AY(t.) — H(t)A©(*) ||, is greater than some threshold R. R is set to meet a 
continuity requirement - that is, not to exceed the allowed number of false alarms of system 
malfunction caused by regular measurement noise. For a given acceptable error radius a, 
we can only guarantee that RAIM will alert us to position errors || r u (t) — r u (t) \\> a 
using the threshold R for given satellite geometries. Fig. 15 addresses the availability of 
such geometries with respect to latitude to alert for radial errors of 1.1m while allowing a 
continuity risk of 2 x 10" 6 per 15 sec. These results were for 122.17 deg West. We assume 
that GPS is augmented only with the Globalstar Constellation. A conservative phase 
noise variance of 1.4 cm, and phase reading rate of 5 Hz are assumed. 

VIILC. A Joint GPS /Globalstar Reference Transceiver 

Each Globalstar satellite downlink contains 13 1.23 MHz bands, spanning 2483.5 MHz 
to 2500 MHz in the S-band. Each band supports 128 CDMA channels, one of which 
harbors a pilot signal, the modulation of which can be described by equ. (1). In this 
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context, D(t) refers to an outer PN sequence of length 288 which is chipped at l.2Kcps. 
Cj(t) is created from the sum of 2 inner PN sequences, each of length 2 10 . which is then 
filtered to a 1.23MHz bandwidth. Cq(£) is a different code, created in a similar manner 
to Ci(t). Fig. 16 illustrates an architecture for the joint GPS and Globalstar reference 
receiver 3, together with a transmitter subsection 90. The rf filter 23 has a center 
frequency of f G — 2492MHz and a bandwidth of roughly 75MHz. The f Grf 24 and 
fdf 28 mixing frequencies can be generated using an integer-N synthesizer 89 with a 
dual modulus prescaler (see Lee [26]). The if mixing scheme 27 adheres to that of Fig. 
4b. The signal is filtered 36 to a pre- correlation two-sided bandwidth, B c of less than 
% 2.5MHz. In-phase and quadrature sampling 37 of the Globalstar signal occurs at 
f s « 20MHz, the clock speed 99 of the Globalstar tracking assembly 88 circuitry. In 
the preferred embodiment, the LRU 8 is implemented using a ground-based transmitter 
90 which employs a VHF rf carrier frequency f t . The f t carrier signal 96 is generated 
by the frequency synthesizer 89. The data modulation section 91 is controlled by the 
microprocessor 56 via the data/control bus 93. Data is modulated onto the carrier using 
PSK modulation at a rate of 9600 baud. The signal is then amplified 94 and transmitted 
via a VHF antenna 95. 

VIII. D. Primary Error Sources for Precise Navigation 

This section gauges those additional sources of error which cannot be directly estimated 
by the data reduction algorithm of section VI. A. We provide only enough detail to roughly 
calibrate the sources of error. 

VIII. D.l Receiver Phase- Tracking Errors 

We will gauge this error using equ. (9), since the filtering of the Globalstar inner PN 
sequences does not have a substantial effect on the phase tracking performance. It is 
prudent, for most implementations of the invention, to select a narrow phase-locked loop 
bandwidth Bl « 10 Hz. This Bl enables a second-order phase-locked loop with damping 
ratio C 5=3 0.7 to track phase acceleration of 100.27rrad/sec 2 with an error < O.lrad. 
For a receiver noise figure of roughly 3dB> the nominal Globalstar transmission achieves 
^» = 37.5 dB — Hz. At Bl = IQHz, we expect a 1 — a phase error due to thermal noise 
of roughly 0.\2rad. 

VIII. D.2 Ionospheric Errors 

The distance-equivalent group delay due to the ionosphere can be as large 20m. However, 
if a differential carrier-phase measurement is taken as in equ. (20), and the user and 
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reference are within lO&m, the resultant error is governed by local irregularities in the 
ionospheric structure, which delay the signal to the user and reference station by different 
amounts. For S-band transmissions, and a user- reference separation of d km, we estimate 
the resultant phase error as a normally distributed random process of zero mean and 
variance A 2 4.4 x 10" 4 d rad 2 , where A is the wavelength in cm. This leads to 1 — a phase 
deviations on Globalstar signals of 0.25rad and 0.57 rad for distances of Ifcm and bkm 
respectively. The corresponding GPS deviations are 0.44 rad and 0.99 rad respectively. 

VIII.D.3 Tropospheric Errors 

Without any form of differential correction, the delays caused by the troposphere at a 
satellite elevation of 10 dcg are roughly 14m. With differential measurement, and a baseline 
separation of d km < 10 km the remaining tropospheric delay can be roughly modeled as 
a normal distribution ~ N(0 cm,(0.1d) 2 cm 2 ). 

VIII.D.4 Ephemeris Errors 

For the Globalstar satellites, we expect the ephemeris disturbances discussed in section 
(VI. A) to be bounded by 1.5 cm and 7.2 cm for d of 1 km and 5 km respectively. The 
ephemerides of the GPS satellites are known to within roughly 10 m rms; the resulting 
ephemeris errors are bounded by 0.05 cm and 0.25 cm respectively. 

VIII. E. Expected performance of a System Using only the Globalstar Constellation 

To illustrate the performance of the invention, we discuss a Monte-Carlo simulation which 
indicates the behavior of a system augmenting GPS with the Globalstar Satellites. Sepa- 
rate simulations are conducted for separations between user and reference station of lArm 
and 5Jfcm. The simulations are all conducted assuming the user is in Palo Alto, California, 
and is capable of seeing all satellites above 10 dcg elevation. For this study, (see Rabinowitz 
[27]), conservative disturbances were assumed due to ionospheric and tropospheric delay, 
thermal noise, ephemeris errors, and oscillator instabilities. For each simulation, the as- 
sumed time of the experiment was varied, sequentially sampling the interval of 12 hours, 
roughly one period of the Navstar satellite orbits. For each simulation, the user receiver 
is ascribed velocity and displacement from the reference receiver in a random direction. 
It was assumed that the user's motion was relatively slow, so that the user-reference 
separation was roughly constant over the course of tracking. 

Fig. 17 displays the 1 — a deviation of radial position errors as a function of tracking 
time for a mobile user. The parameter estimates were found using equ. (59). Each point 
corresponds to the mean error deviation averaged over 200 simulations. Fig. IS displays 
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the evolution of the lower bound on the probability of selecting the correct set of integers 
for the Navstar satellites in view over the tracking period. Each point in this figure 
represents the worst case for 200 simulations. This lower probability bound is calculated 
according to the technique discussed in section (VI.C.4). Once the integers are correctly 
identified, a user may rely completely on GPS measurements where ephemeris errors are 
negligible. Since the measurement errors for each satellite are roughly distributed as 
~ Af(0,a£ s ), the positioning error deviation for a static user decreases roughly as oc 
where N is the number of measurements taken on each satellite. 
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IX. Claims 

What is claimed is: 

1. In a satellite-based navigation system, a method for estimating a precise position of a 
user device, the method comprising: 

transmitting carrier signals from a set of satellites, wherein the set of 
satellites includes a set of LEO satellites; 

sampling at a reference station the carrier signals to obtain reference carrier 
phase information comprising geometrically diverse reference carrier phase information 
from the set of LEO satellites; 

sampling at a user device the carrier signals to obtain user carrier phase 
information comprising geometrically diverse user carrier phase information from the set 
of LEO satellites; 

calculating the precise position of the user device based on the reference 
carrier phase information and the user carrier phase information, wherein the calculation 
uses the geometrically diverse reference and user carrier information from the set of LEO 
satellites to quickly resolve parameters related to the integer cycle ambiguities in the 
reference and user carrier phase information. 

2. The method further comprising: 

receiving code signals from a set of navigation satellites; 

measuring at the reference station code signals to obtain jeference code 
phase information; 

measuring at the user device the code signals to obtain user code phase 

information; 

estimating user and reference clock biases from the user and reference code 
phase information; 

correcting for clock offsets using the estimated user and reference clock 

biases. 

3. The method further comprising: 

estimating an approximate user position using code phase signals received 
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from a set of navigational satellites. 

4. The method further comprising: 

communicating differential code phase correction data and the reference 
carrier phase information from the reference station to the user device. 

5. The method further comprising: 

communicating LEO satellite ephemeris data to the user device directly 
from the reference station, or using a satellite data link. 

6. The method wherein the step of calculating the precise position of the user device 
comprises predicting present reference carrier phase information based on past reference 
carrier phase information. 

7. The method wherein the step of calculating the precise position of the user device com- 
prises compensating for frequency dependent phase delay differences between navigation 
carrier signals and LEO carrier signals in user and reference receiver circuits. 

8. The method wherein the step of reading the carrier phase information at the user device 
comprises reading navigation carrier information and LEO carrier information within a 
predetermined time interval selected in dependence upon an expected motion of the user 
receiver and the motion of the LEO signal sources. 

9. The method wherein the calculating step comprises correlating carrier phase informa- 
tion from a first LEO beam with carrier phase information from a second LEO beam. 

10. The method wherein the calculating step comprises calibrating LEO oscillator insta- 
bilities using navigation satellite information. 

11. The method wherein the calculating step comprises compensating for phase distur- 
bances resulting from a bent pipe LEO communication architecture. 

12. The method further comprising the step of 

monitoring the integrity of the calculating step. 

13. A satellite-based navigation system comprising: 

a set of satellites transmitting carrier signals, wherein the set of satellites 
includes a set of LEO satellites; 

a reference station tracking the carrier signals to obtain reference carrier 
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phase information comprising geometrically diverse reference carrier phase information 
from the set of LEO satellites; 

a user device; and 

a communications link between the reference station and the user device; 
wherein the user device comprises 

a receiver tracking the carrier signals to obtain user carrier phase informa- 
tion comprising geometrically diverse user carrier phase information from the set of LEO 
satellites; and 

a microprocessor calculating a precise position of the user device based on 
the reference carrier phase information and the user carrier phase information, wherein 
the calculation uses the geometrically diverse reference and user carrier information from 
the LEO satellites to quickly resolve parameters related to integer cycle ambiguities in 
the reference and user carrier phase information. 

14. The system wherein the set of satellites comprise a set of navigational satellites, and 
wherein the communications link conveys differential code phase correction data and the 
•reference carrier phase information from the reference station to the user device. 

15. The system wherein the communications link conveys LEO satellite ephemeris data 
to the user device from the reference station. 

16. A user device for providing satellite- based navigation, the device comprising: 

at least one antenna for coupling to carrier signals transmitted from a set 
of satellites, wherein the set of satellites includes a set of LEO satellites; 

a first receiver for tracking the carrier signals to obtain user carrier phase 
information comprising geometrically diverse user carrier phase information from the LEO 
satellites; 

a second receiver, not necessarily distinct from the first receiver, for re- 
ceiving reference carrier phase information transmitted from a reference station, wherein 
the reference carrier phase information comprises geometrically diverse reference carrier 
phase information from the set of LEO satellites; 

a microprocessor for calculating the precise position of the user device based 
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on the reference carrier phase information and the user carrier phase information, wherein 
the calculation uses the geometrically diverse reference and user carrier information from 
the set of LEO satellites to quickly resolve parameters related to integer cycle ambiguities 
in the reference and user carrier phase information. 

17. The device wherein: 

the set of satellites comprises navigation satellites; 

the first receiver measures navigation code signals to obtain user code phase 

information; 

the second receiver receives reference code phase information transmitted 
from the reference station; and 

the microprocessor estimates user and reference clock biases from the user 
and reference code phase information, and uses the estimated clock biases to correct for 
clock offset errors. 

18. The device wherein the first receiver reads navigation carrier phase information and 
LEO carrier phase information at times separated by no more than a predetermined time 
interval selected in dependence upon expected movement of the device and the movement 
of the LEO satellites. 

19. A user device for providing satellite-based navigation, the device comprising: 

at least one antenna for coupling to signals transmitted from a set of 
satellites, wherein the set of satellites includes a set of navigation satellites and a set of 
LEO satellites; 

a receiver for tracking the signals to obtain code phase information and 
carrier phase information comprising geometrically diverse carrier phase information from 
the LEO satellites; and 

a microprocessor for calculating the precise position of the user device based 
on the code phase information and the carrier phase information, wherein the calculation 
uses the geometrically diverse user carrier information from the LEO satellites to resolve 
parameters related to integer cycle ambiguities in the carrier phase information from the 
navigation satellites. 
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